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SUMMARY 

To be feasible for computationally intensive applications such as parametric studies, optimization and 
control design, large-scale finite element analysis requires model order reduction. This is particularly 
true in nonlinear settings that tend to dramatically increase computational complexity. Although 
significant progress has been achieved in the development of computational approaches for the 
reduction of nonlinear computational mechanics models, addressing the issue of contact remains 
a major hurdle. To this effect, this paper introduces a projection-based model reduction approach 
for both static and dynamic contact problems. It features the application of a non-negative matrix 
factorization scheme to the construction of a positive reduced-order basis for the contact forces, and 
a greedy sampling algorithm coupled with an error indicator for achieving robustness with respect to 
model parameter variations. The proposed approach is successfully demonstrated for the reduction 
of several two-dimensional, simple, but representative contact and self contact computational models. 
Copyright © 0000 John Wiley & Sons, Ltd. 
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1. INTRODUCTION 

The nonlinear finite element (FE) analysis (FEA) of large-scale systems often requires 
prohibitively large computational resources. These are typically prescribed by the fine 
discretization of the computational domain that leads to a large number of degrees of freedom 
(dofs) in the system, as well as, in the case of dynamic analysis, the potentially large number 
of time-steps needed to accurately describe the evolution of the system. 

Projection-based model order reduction (MOR) techniques alleviate the first issue by 
restricting the solution space to a smaller subspace, thereby reducing the number of dofs. 
While many approaches have already been developed for the efficient reduction of linear 
computational models 13 II H H, three main strategies have been explored so far for 
efficiently reducing nonlinear computational models. The first one is based on linearization 
techniques 5, 6, The second one is based on the notion of pre-computations mm but is 
limited to polynomial nonlinearities. The third strategy relies on the concept of hyper-reduction 
■ that is, the approximation of the reduced operators underlying a nonlinear reduced- 
order model (ROM) by a scalable numerical technique based on a reduced computational 
domain Eunannnaiimiiaiia- While several hyper-reduction techniques have been proposed 
in the literature, two of them have been designed specifically for the nonlinear FEA of structural 
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and solid mechanics problems: The a priori hyper-reduction method [9], and the Energy 
Conserving Sampling and Weighting method [10, T4j [16]. Nevertheless, contact problems 
remain a major hurdle for nonlinear model reduction in the context of structural analysis. This 
is because, among other things, contact problems are characterized by inequality constraints 
that complicate the reduction process. 

Most if not all (nonlinear) computational contact methods proceed in two steps. The first 
one focuses on contact detection — that is, the identification of nodes, edges, and/or faces 
of the computational model that should be in contact. The second step focuses on contact 
enforcement — that is, the satisfaction of the contact constraints defined by physical laws 
such as non-penetration and frictional behavior. In practice, the aforementioned constraints are 
enforced using one of three popular approaches: The penalty method, the Lagrange multiplier 
method, or the augmented Lagrange multiplier method mm- Attention is focused in this 
paper on the case of the Lagrange multipler method (or its augmented version). 

For large-scale nonlinear dynamical systems, proper orthogonal decomposition (POD) [IjJj 
is the method of choice for generating the reduced-order basis (ROB) needed for constructing 
a ROM. It proceeds by collecting solution snapshots during a training procedure, then 
compressing them using the singular value decomposition (SVD) method. The resulting ROB 
minimizes the projection error of the snapshots. However, when the Lagrange multiplier 
method is chosen for enforcing contact constraints, reducing the contact forces — or 
equivalently, the Lagrange multipliers — requires special care because the reduced multipliers 
must have a positive sign. A ROB for these dual variables based on SVD would not enforce this 
positivity requirement a priori. For this reason, it is proposed in this paper to reduce positive 
quantities such as contact forces or Lagrange multipliers using a positive counterpart of the 
SVD method known as the non-negative matrix factorization (NNMF) method. Introduced 
first in the context of image compression [2D], this method builds low-rank positive factors 
that approximate a given non-negative matrix. Here, it is shown that the proposed usage 
of NNMF results in the construction of a ROB that can accurately represent the Lagrange 
multipliers and lead to the effective reduction of contact computational models. 

In [21, [22], the authors have addressed similar issues by constructing a ROB for the 
Lagrange multipliers using a positive linear combination of pre-computed snapshots of these 
dual variables. For time-dependent problems, this approach can rapidly become impractical as 
it can lead to the construction of a ROB of very large dimension. In this work, NNMF provides 
a natural procedure for optimally compressing a potentially large number of snapshots of the 
dual variables and constructing a small dimensional ROB for approximating them. 

A more generic model reduction issue is the robustness of a ROM with respect to variations 
of the model parameters. Indeed, a ROM is truly useful when it can be used as a surrogate of 
the underlying high-dimensional model (HDM) for parameter values that may be different from 
those sampled for the purpose of constructing a ROB. Contact problems can be particularly 
sensitive to parameter variations, for example, when the contact areas are very sensitive to such 
variations. A popular approach for constructing a ROM that is valid in a large region of the 
model parameter domain is to couple a greedy approach with one or several a posteriori error 
estimators mm mm in order to effectively sample the parameter domain for computing 
solution snapshots and constructing a ROB. Such an approach constructs increasingly accurate 
ROMs by detecting locations in the model parameter domain where the errors associated with 
the ROM are the largest. The associated HDM is subsequently reconstructed at the identified 
worst-error parameter values and solution snapshots are computed and stored. Then, the ROB 
and therefore ROM — is updated based on these additional snapshots, thereby reducing 
drastically the error(s) for the newly sampled parameter values. This procedure terminates 
whenever the estimated ROM error(s) is (are) below a specified tolerance throughout the 
model parameter domain of interest. It leads to a ROM that remains accurate away from 
the training configurations. Therefore in this work, a greedy approach is also developed to 
construct both primal and dual ROBs that are robust in a given model parameter domain. 
Specifically, the greedy approach developed in this paper relies on the definition of an error 
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indicator for the contact problem, and successive updates of both primal and dual ROBs 
computed using SVD and NNMF, respectively, are performed. 

The remainder of this paper is organized as follows. The notation adopted in this paper 
is presented in Section [2] The considered family of contact problems is derived in Section [3] 
The proposed model reduction procedure for contact problems is described in Section [4] This 
procedure is applied in Section [5] to the reduction of three different contact problems. Finally, 
conclusions are offered in Section [g] 


2. NOTATION 

Throughout this paper, matrices are denoted by bold capitals (ex. A), vectors by bold lower 
cases (ex. a), and subscripts identify rows and columns (ex. A^j is the entry of A located in 
the i-th row and j-th column of this matrix). 

A + denotes the Moore-Penrose pseudo-inverse of the matrix A. 

In identifies the identity matrix of size N and 0 identifies a matrix of zeros. In identifies 
the vector of dimension N whose entries are all ones. 

For two matrices A and B of equal dimension M x N , the Hadamard product A © B is 
the matrix of the same dimension whose entries are given by 

(A © B) it j = Aij ■ Bij (1) 

A discretized variable u £ at time-step n £ N is identified by a superscript as u n £ l w . 
The standard Euclidean norm of a vector x £ and the Frobenius norm of a matrix 
A £ K Mxiv are denoted by ||a;|| 2 and ||A|| F , respectively, and defined as follows 

l l 

( N \ 2 / M N \ 2 

, I|A|| f = ■ ( 2 ) 

i -1 / \1=1 i=l / 

Finally, the negative part of a real number x is defined as [x]_ := min(;r,0) and that of a 
vector x £ R w is defined as [cc]_ = [[xj]_]. 

3. MODEL CONTACT PROBLEM 

The main issue addressed in this paper is that of how to reduce the dual Lagrange multipliers 
introduced in a solution process for enforcing the inequality constraints governing a static or 
dynamic contact problem. For this reason, and for the sake of clarity, attention is focused 
here on a model contact problem where the individual bodies in contact (for example, see 
Fig. 0 have linear material and kinematic behaviors. The reduction of the primal displacement 
solution in the presence of material nonlinearities and/or large displacements and rotations 
raises independent issues that have already been addressed elsewhere in the literature, for 
example, in mmmm- Furthermore, for simplicity and without any loss of generality, 
only the case of a frictionless, adhesive-free normal contact is considered, and all body problems 
are assumed to be undamped and semi-discretized on uniform matching meshes. 

Using the modeling assumptions stated above, the FE semi-discretization of a dynamic 
TVo-body contact problem can be written in matrix form as 


Mil + Ku = f + B t A 
Bu — c > 0 

•u(O) = u 0 

u( 0 ) = u 0 


(3) 
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where a dot designates a time derivative, the first semi-discrete equation expresses the dynamic 
equilibrium of Nq given (flexible) bodies, the inequality constraint derives from the semi¬ 
discretization of the Hertz-Signorini-Moreau contact conditions that coincide with the Karush- 
Kuhn-Tucker (KKT) complementary conditions in the theory of optimization [13], and the 
last two equations result from the semi-discretization of the initial conditions. In Eq. ([3| 
above, M £ R NxN and K £ R NxN , where N denotes the total number of primal dofs, are 
constant symmetric positive definite or semi-definite block diagonal mass and stiffness matrices, 
u = u(t) £ is a semi-discrete displacement vector and t denotes time, / = f(t) £ l w is 
a semi-discrete force vector, B £ R NxXN , where N\ denotes the total number of potentially 
active contact inequalities, is a signed boolean matrix which extracts from u the pairs of dof 
governed by a contact condition, c is the vector of initial clearances, and A is the vector of 
semi-discrete Lagrange multipliers. 

Specifically, M, K , u, B , A, and f can be written as 



~M 1 

0 

0 



0 

0 


Mi 


0 

m 2 ■ 

0 


0 

K> 

0 


«2 

M = 

0 


• M Nn_ 

, K = 

0 


1 

... 

, u = 

u N n 





Ai 


' fl ' 




^2 


f 2 

B = [B ± B 2 

' B N n ] , 

A = 

_^N n _ 

, f = 

fN n _ 


(4) 


where the subscript i = 1,2, • • ■ , Nq designates the body f 

Again, for simplicity and without any loss of generality, an implicit time-discretization is 
assumed in the dynamic case so that solving both static and dynamic contact problems can 
be formulated as 


( u n , A") = argmin -v T Av — v T b n — fi T {Bv — c) (5) 

uei». M eR+ JVx 2 

where T designates the transpose operation, the superscript n designates the n-th time-step 
t n , to = 0 < fi < • • ■ < t n < ■ ■ ■ < tN t = T is a discretization of [0, T] into N t + 1 time-points, 
A £ R NxN is the block diagonal matrix containing for each body its mass and stiffness matrices 
scaled according to the chosen time-integration algorithm and selected time-stepping strategy, 
b £ contains for each body the right-hand side vector arising from the implicit time- 
discretization of the dynamic semi-discrete equations of equilibrium governing this body, and 
A > 0. In the static case, the superscript n is dropped, A = K, and b = f. 

4. MODEL REDUCTION 


4-1. Galerkin projection 

For the model contact problem described above, the standard Galerkin projection-based MOR 
method is appropriate. In this method, the primal and dual components of the solution 
are approximated in two reduced subspaces represented here by two pre-computed ROBs 
U £ R Nxp and U\ £ R N * xpx , respectively. This can be written as 

«"( 7 ) * Uu n r { 7), A" (7) « U x X n r (7) (6) 

where u r E and A r E MP X are the generalized coordinates of the reduced displacement and 
Lagrange multiplier solutions, respectively, and 7 E V C M m is a vector of m parameters of 
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the contact problem of interest. Inserting the above two approximations in the saddle point 
problem ([5| gives 


(it", A") = argmin -v r T A r v r — v r T b r — n r T (B r v r — c r ) (j'j 

u r £H>, n r es.+ FX ^ 

where 

A r := U T AU G M. pxp , K ■= U T b n G R p , B r := U^BU G R pxXp , and c r := f/Jc G R Px . 
To ensure the non-penetration condition in the contact ROM, the reduced vector of Lagrange 


for satisfying this requirement is to construct a non-negative ROB U\ > 0. Indeed, since the 
solution of the contact problem formulated using the contact ROM delivers a positive vector 
of generalized Lagrange multiplier coordinates A, >0, U\ > 0 guarantees in this case that 

U\X r > 0. 

Algorithm [l] below outlines the Galerkin projection-based MOR method proposed in 
this paper for solving the contact problem ([5|. To this effect, note that in general, it is 
feasible to compute the reduced vector fo" online. Specifically, this can be achieved by 
pre-computing some relevant small-size quantities offline. For example, consider the case 
where the prescribed, time-dependent force vector can be decomposed as f(t) = Lg(t), 
where L G l w describes the time-invariant spatial distribution of / and j(t) G t describes 
its temporal evolution. If time-discretization is performed using the midpoint rule, 
b n r = M r < + (At/2)M r ti" + (At 2 /4)U T Lg(t n+1/2 ), where At is the computational time- 
step. In this case, Iff can be efficiently computed at each time-step by pre-computing once for 
all the quantity (At 2 /A)U T L. 


multipliers must be non-negative — that is, £/\A" > 0. The approach proposed in Section 4.3 


Algorithm 1: Online solution of the contact problem i> using a Galerkin projection-based 
contact ROM _ 

input : Reduced quantities A r , B r , c r , uff and iff 
output: Generalized coordinates 1 1 , {A"}^ r i 1 


1 

2 
3 


for n = 1,2,..., N t do 

Construct the reduced vector 6"; 

Solve the reduced saddle point problem |7j) 


(tt",A")= argmin 

Ur-eRp, /i r eR+ PX 



A-f.Vj. - Vy. by. 


/x r T (B r v r 


Cr) 


4 end 


REMARKS. 

• In general, the primal ROB U enjoys an orthogonality property, and the initial conditions 
are specified for the high-dimensional fields u° and iff which have physical meanings. 
Hence, these initial conditions can be converted into initial conditions for the generalized 
coordinates iff and iff using ([6| and the orthogonality property of U. 

• Both dual unknowns A and A r are auxiliary variables. They do not necessarily need 
special initializations. 

4-2. Construction of an optimal primal reduced-order basis 

For the model contact problem described in Section [3j the adoption of global ROBs for 
both primal and displacement components of the solution is appropriate. More complex 
contact problems featuring material and/or geometric nonlinearities within the solid bodies 
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call however for adaptive primal and dual ROBs. Such ROBs can be constructed, for example, 
using the concept of locality introduced in [13] which does not necessarily refer to space or 
time, but to the region of the manifold where the nonlinear solution lies. 

In either case, a primal (dual) ROB can be constructed from the compression of primal 
(dual) solution snapshots — that is, primal (dual) components of solutions of problem ([5| for 
different time-instances tj and different instances 7 S of the parameter vector 7 . Specifically, 
for each sampled parameter vector 7 ^, s = 1, • • • ,N S , the computed primal and dual snapshots 
are gathered in matrices X s and Xj[, respectively, with Xf^ := uj(y s ) and X%. . := Aj(-y s ), 
j = 0,---,N t . 

In general, the primal ROB is not subject to any particular constraint. Therefore, it can 
be constructed by POD m via the SVD decomposition of the global primal snapshot matrix 
X := [X 1 ,..., X Ws ]. This corresponds to solving the optimization problem 


minimize IIX — UV\\ 2 F 

to compute the low-rank approximation of X 

X ss UV 


( 8 ) 

(9) 


Hence, the ROB U is constituted of the first p left singular vectors of the snapshot matrix 
X and V = T,W t , where S is the diagonal matrix of the first p singular values of X, and W 
is the matrix of its first p right singular vectors. 


4-3. Construction of an optimal dual reduced-order basis 

As emphasized in Section [4T| it is essential to preserve the positivity of the contact constraints 
after reduction. For this purpose, it is proposed here to construct a positive dual ROB U\ using 
NNMF [20]. This corresponding to solving the optimization problem 

minimize \\X x - U\V\\\ 2 F 

u x eR w * x , v A eR^x x N t 

subject to U\ > 0 (I®) 

Va>0 

where X A := [X{,...,X^] is the global dual snapshot matrix. The NNMF algorithm leads 
to the low-rank approximation of the dual global snapshot matrix by two positive factors 


Unlike problem 


*A « U\V\ (11) 

problem (10) does not have a closed form solution. Consequently, 
this problem is usually solved using an iterative method that typically converges to a local 
minimum. Examples of such methods are the original multiplicative updating rule [2D], 
the alternating non-negativity least-squares method EH, and block coordinate descent 
algorithms [28] . 


4-4- Snapshot selection 

When many snaphsots are collected for the purpose of constructing primal and dual ROBs 
with a potential for accurate approximations of the displacement and Lagrange multiplier 
fields, respectively, the SVD and NNMF of the global snapshot matrices X and X\ become 
computationally intensive. When N s different instances of the parameter vector 7 are sampled, 
both aforementioned matrices have N s x (7V t + 1) columns. In this case, the following remarks 
are noteworthy: 

1 . It is not necessary to store the dual snapshot solutions at those time-steps where there is 
no contact, as these snapshots are zero. Hence, the dimension of the matrix X\ can be 
reduced to the number of time-steps at which contact is established, without modifying 
the accuracy of the resulting ROM. 
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2. The dimensions of the global snapshot matrices can be further reduced by down-sampling 
the HDM in time. This may be necessary when a very large number of time-steps N t 
is computed, for example, when the solution snapshots are obtained using explicit time¬ 
stepping in a relatively large time-window [0,7”]. The temporal down-sampling of the 
solution snapshots may affect however the accuracy of the resulting ROM as it implies 
fewer information for training this ROM. 

In the remainder of this paper, temporally down-sampled sets of primal and dual solution 
snapshots are denoted by {u n } ne 5 and {A n }„ e s A , respectively, where S C {0, ■ • ■ ,N t } and 
S\ C {0, • • • , N t }. In Section [73] a preliminary study of the effect of temporal down-sampling 
of the solution snapshots on ROM accuracy is performed. 

4-5. Construction of parametrically robust ROBs 

In summary, primal and dual ROBs can be constructed by compressing primal 
and dual components of solution snapshots computed for some parameter instances 
7 S GT>, s = 1, • • • ,N S . To this effect, it is first noted that an a priori sampling of the 
parameter space may miss certain regions of T> where the ROM will be inaccurate. This 
underscores the importance of sampling 7 at specific instances 7 S that enable the construction 
of a parametrically robust ROM — that is, a ROM that is accurate in the entire parameter 
space V. 

Finding the best samples in T> is however a combinatorial problem whose solution is 
often intractable. For this reason, economical greedy strategies have been developed for this 
purpose j23j EH EH ■ Such strategies proceed iteratively by identifiying the parameter samples 
for which the error associated with the current ROM is the largest, then sampling the HDM at 
these parameter instances, and finally updating the ROM using the additional HDM solution 
snapshots. 

Finding parameter instances 7 S that maximize the error of the current ROM can 
be performed by solving directly an error maximization problem using a gradient-based 
optimization algorithm [24], or a global optimization approach and a surrogate model [25] . 
Alternatively, a basic greedy procedure [23] can be designed for this purpose as follows. 
Given an a priori set of N c candidate parameter instances { 7 * 1 \--- , 7 ^°)} C 2?, where 
the superscript and pair of parentheses emphasize here the candidate aspect of a parameter 
instance 7 W and distinguish it from an effectively sampled parameter instance 7 S , choose the 
elements of this set which maximize the norm of the error 


|e(7)lII := 


II UK 


(7)-« n (7)ll^ 


\n—0 


( 12 ) 


between the HDM solution u"( 7 ) and the associated ROM solution Uu™(j). In practice 
however, the set of HDM solutions {w"(7)}^lo is unknown. Therefore, the above error is 
conveniently replaced with an error indicator (that is preferrably economical). Here, this 
indicator is based on the following contact conditions: 


(B+(Au n 


Bu n - c > 0 

b n )) © ( Bu n - c) = 0 
B+{Au n -b n ) >0 


(non-penetration) 
(complementary slackness) 
(contact force positivity) 


(13) 


Specifically, given a set of ROM solutions {v,™ (7)}^1 0 , th e error indicator proposed in this 
paper is 


1(~r,a lt a 2 ,a 2 ) ■= £ (a^«( 7 )) 2 + a 2 ||rj( 7 )||l + (14) 

nej 
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where 


r i (7) = B{l)Uv%{~f)-c(~i) 

r 2 ( 7 ) = ( B+ {{'Y)( A (,l) u K('r) - b n h))) © (B(7)i7M^(7) - 0(7)) (15) 

^ 3 ( 7 ) = s+ (7)(^(7)^<(7) - &"(7)) 


.7 defines a subset of time-steps at which the proposed error indicator is evaluated, and 
(j>{v) := 11[tj]_ 11 2 - The coefficients a,, z = 1, • - ■ ,3 are adjustable weights that can be used 
to emphasize the relative importance of each contact condition. 

Because 1 ( 7 , aq, a 2 , < 33 ) characterizes the violation of the contact conditions, it is an 
indicator for the error associated with the ROM solution Uu r { 7 ). The computational 
complexity of this error indicator is reasonable in the sense that its evaluation does not require 
the solution of any system of equations. It requires only multiplications. Of particular interest 
is the case «2 = 0 for which the evaluation of 1 ( 7 , oq, 0, a 3 ) becomes the most economical. In 
fact, the applications discussed in Section[5]suggest that the case aq = 1, ct 2 = «3 = 0 leads to a 
good error indicator. As for the time-sampling of the snapshots, the case J = {0, • • • , N t } leads 
to the most accurate error indicator. Unfortunately, it generates an excessive computational 
burden when N t is very large. 

In any case, the reader is reminded that parameter sampling is an integral part of a training 
procedure that is performed offline. Therefore, the computational complexity of the error 
indicator (141 does not affect the performance of the ROM computations to be performed 
online. 

The training procedure adopted in this work is summarized in Algorithm [2| While it relies 
on the basic greedy approach outlined above, this procedure can be accelerated using the 
techniques described in (SB]. 


5. APPLICATIONS 

The model reduction approach proposed in this paper for contact problems is illustrated here 
with three simple but representative two-dimensional, parameterized, model problems. The 
first one is a static problem of the obstacle type. The two other ones are dynamic contact 
problems. In the first application, the Lagrange multipliers are approximated using a positive 
linear combination of the computed snapshots. Data compression is not performed in this case 
because the number of snapshots computed during the training procedure is small. Hence, 
this first problem is designed to demonstrate in particular the effectiveness of the proposed 
greedy algorithm (Algorithm [ 2 ]) for sampling the parameter space. The second model problem 
considered herein is a dynamic version of the first problem. Its dynamic aspect gives the 
opportunity to precompute a large number of solution snapshots. Hence, it is suitable for 
illustrating the effectiveness of the proposed approach for constructing a ROB for the Lagrange 
multiplier field. The third considered model problem is a dynamic contact problem between two 
parallel Kirchoff plates. It serves the purpose of demonstrating the applicability of the proposed 
model reduction approach to more generic, parameterized, multi-body contact problems. For 
the first two problems, the performance of a constructed ROM is assessed in “predictive mode” 
that is, for a problem configuration different from that used for training the ROM. For the 
third problem, the performance of a ROM is assessed in “reproduction” mode — that is, for 
the same problem configuration as that used for training the ROM. 

In all cases, the relative error of the solution delivered by a ROM is defined in the static 
case as 


relative error (%) := 


Uuri'j) - u(-y)\\l 
ll«(7)lli 


X 100 


(16) 
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Algorithm 2 : Greedy sampling algorithm for dynamic contact problems 
input : Initial sampled parameter instance 71 , set of N c candidate parameter instances 
C = { 7 « ...,7 (iVo )} C T> 7 maximal number of primal and dual basis vectors p 
and p\ 7 respectively, maximum number of greedy iterations N gre edy > 2 , 
convergence tolerance e < 1 

output: Global primal and dual ROBs U and U\, respectively, number of sampled 
parameter instances N s 


1 

2 

3 


for Niter — I; * • * ; A greedy dO 

Compute HDM snapshots {u n (j jv itel .)}neS and {A”( 7 jv iter ,)} raei s A and store them in 
local snapshots matrices X NitBr and X^ iter 

Accumulate global primal and dual snapshot matrices X and X\, respectively 


X = [A 1 , • • • , X NitBr ] , X x =[Xl 


4 

5 

6 

7 

8 
9 

10 

11 


12 

13 

14 

15 

16 

17 

18 
19 end 


Construct a primal ROB U of dimension < p by compressing X using SVD 
Construct a dual ROB U\ of dimension < p\ by compressing X\ using NNMF 
if Niter = Ng r eedy then 

Ng — Ngreedy 

terminate the algorithm 
end 

for j = 1,..., N c do 

Compute ROM solutions {«” ( 7 ^)}^L 0 using Algorithm [T] and the current ROBs 
U and U x 

Compute the a posteriori error indicator X ( 7 ^) 

end 

Find 7 iv iter +i = argmax 7eC J( 7 ) 
if X(j Niter+1 ) > 61(72) then 
N s = Niter 

terminate the algorithm 
end 


where u( 7 ) is the static HDM solution, and in the dynamic case as 


relative error {%)'■= 


N t 


E ||E/u?(7)-u"(7)II! 

n —0 


N t 


x 100 


E ll M ”(7)ll! 

n —0 


where u n { 7 ) is the dynamic HDM solution. 


(17) 


5.1. Static problem of the obstacle type 


The model problem presented here is that of the computation of the equilibrium position 
of a two-dimensional elastic membrane covering the spatial domain (x, y) G [ 0 , 1 ] x [ 0 , 1 ], 
constrained by homogeneous Dirichlet boundary conditions along all its boundaries, subjected 
to a uniform load / = —10, and facing a parameterized obstacle. This model, parametric, static 
contact problem can be described by the inequality-constrained Poisson equation 


V 2 u = f 

u > girt) 


(18) 
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where 


ff ( 7 ) = -1 + 0_ 4e -200((r C - 7l ) 2 + (y-0.5) 2 ) + 72e -355.56((x-0.7) 2 +(?/-0.5) 2 ) ( 19 ) 

describes the parameterized obstacle. The range of interest of the two-dimensional parameter 
domain V = ( 71 , 72 ) is set to [0.3,0.6] x [0.2, 0.6]. 

The elastic membrane is discretized into 200 x 200 finite elements, which generates an HDM 
with 40,000 dofs. For the training procedure, V is initially sampled on a 10 x 10 uniform tensor 
grid to generate a set of N c = 100 candidate parameter instances. Then, Algorithm [2] is applied 
to construct two global primal and dual ROBs U and U\ , respectively, and their associated 
ROM. Because at each iteration of the number of computed snapshots is much smaller than 
N c = 100, these snapshots are not compressed. Instead, they are directly used to gradually 
construct U and U\. In other words, for this problem, p and p\ are evolved in Algorithm 
[2] as p = k = N iter . Additional ROMs are also built by randomly sampling (a priori ) T> , for 
the purpose of comparing their performance to that of the ROM delivered by Algorithm 
[ 2 ] Specifically, 20 instances ( 71 , 72 ) € V are generated using the Latin Hypercube Sampling 
(LHS) method, and primal and dual ROBs of dimension p = p\ = 20 are constructed using the 
snapshots computed at these sampled parameter values. Because of the randomness associated 
with the LHS samples, this construction process is repeated 50 times. 

Figure[3]reports on the convergence of the greedy sampling algorithm for this model problem. 
The reader can observe that at least in this case, the relatively economical error indicator 
1 ( 7 , 1 , 0 , 0 ) performs well and better than all other considered canonical configurations of this 
indicator. The reader can also observe that all intermediate ROMs constructed using the greedy 
sampling algorithm outperform the ROMs constructed using the a priori random sampling 
of T>. For example, an intermediate ROM constructed with p = p\ ss 10 using Algorithm 
[ 2 ] is shown in Figure [3] to deliver the same accuracy as a twice as large [p = p\ = 20) 
ROM constructed using random sampling of the parameter space. Furthermore, the ROM 
characterized by p = p\ = 20 and delivered by the same greedy sampling algorithm equipped 
with a = {1, 0,0} is found to be an order of magnitude more accurate than a ROM of equivalent 
dimension constructed using a random sampling of T>. 

Figure [4] showcases the performance of two ROMs obtained after 20 iterations of Algorit hm [2| 
Specifically, it compares for two different parameter combinations (7 = (0.6,0.6) in Figure [4(a)] 
and 7 = (0.330,0.377) in Figure |4(b)| ) one-dimensional slices of the two-dimensional HDM 
and ROM solutions. Note that these two parameter instances are chosen here for performance 
assessment for two reasons: (a) neither of them is part of the set of N c candidate parameter 
instances inputted to Algorithm [2] and therefore neither of them is part of the training of the 
constructed ROM, and (b) 7 = (0.330, 0.377) is determined a posteriori to be the parameter 
instance for which the constructed ROM has the largest error. In both cases, the computed 
HDM and ROM solutions are almost indistinguishable, thereby demonstrating the accuracy of 
the proposed approach for constructing a contact ROM. This accuracy is confirmed in Figure[5] 
which compares the full, two-dimensional HDM and ROM solutions of the considered model, 
parametric, static contact problem for 7 = (0.330,0.377). 


5.2. Dynamic contact problem of the obstacle type 

Next, a dynamic version of the otherwise same model contact problem described in Section [5.1| 
is considered here. It is described by the inequality-constrained initial boundary value problem 


d 2 u 2 

u+f 

u > gin) 


u(x, y : 0) = 0 
|( x , 9 : 0)=0 


( 20 ) 
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where t £ [ 0 , 2 ]. 

The same semi-discrete HDM described in Section |5.1| is adopted for this problem. It is 
discretized in time using the implicit second-order Backward Differentiation Formula (BDF) 
scheme and the constant time-step At = 0.005. The parameter domain T> = (71,72) is set 
again to [0.3,0.6] x [0.2,0.6] and initially sampled on the same 10 x 10 uniform tensor grid 
as before, in order to generate a set of N c = 100 candidate parameter instances. These are 
inputted to Algorithm [2] for performing the training procedure. At each A^ ter -th iteration 
of Algorithm [2] the solution snapshots are collected at each time-step. Thus, the number 
of columns of each of the global primal and dual snapshot matrices grows in this case as 
(2/0.005) x N^er = 400 Ni ter . The data compression strategy is chosen so that the size of each 


98 

of the primal and dual ROBs grows as p = p\ = [2 + (N iter — 1)—J. 

Figure [ 6 ] reports on the convergence of Algorithm [2] for this problem. It reveals that in this 
case, all considered strategies for a = {07,0:2,0:3} perform equally well. After 10 iterations, 
each of these strategies leads to a ROM whose maximum amplitude error is more than 100 
times smaller than that of its counterpart computed at the first iteration. 

Figure [7] and Figu re [ 8 | display the time-histories at three different points in space of the 
solutions of problem (201 computed for two different instances of 5 ( 7 ) using in each case: (a) 
the HDM, (b) the ROM obtained after 10 iterations of Algorithm |2j and (c) a variant of this 
ROM where SVD is used instead of NNMF to compress the dual snapshots collected within 
Algorithm[2] Again, note that both parameter instances 7 = (0.6,0.6) and 7 = (0.300,0.288) 
are selected here for performance assessment because neither of them is part of the training 
of the constructed ROM, and because 7 = (0.300, 0.288) is determined a posteriori to be the 
parameter instance for which the constructed ROM has the largest error. The reader can 
observe that as expected, the SVD-based ROMs deliver a poor performance as they do not 
satisfy the positivity condition of the Lagrange multipliers. On the other hand, the NNMF- 
based ROMs deliver a solid performance. The solutions they deliver track well the HDM 
solutions. This solid performance of the proposed approach for constructing a contact ROM 
is confirmed in Figure [9] which focuses on the entire spatial domain. 


5.3. Two-body dynamic contact problem 

Finally, the two-body dynamic contact problem graphically depicted in Figure[l0]is considered 
here. Each of the two bodies, fR and f2 2 , is a homogeneous, isotropic square plate of edge size 
L = 1 m and thickness h = 1 mm. It is modeled as a linearly elastic Kirchhoff-Love plate, and 
therefore governed by the partial differential equation 

B 2 n 

ph— + D\7 2 \7 2 u-f = 0 (21) 

where p denotes the density, u(x,y : t) denotes the transverse displacement field, f{x,y : t) 
denotes a distributed external force per unit area (pressure), 


D = 


Eh 3 

12(1 - v 2 ) 


( 22 ) 


and E and v denote Young’s modulus and Poisson’s ratio, respectively. The two plates are 
assumed to be made of the same material characterized by p = 7800 kg/m 3 , E = 200 GPa, and 
v = 0.3. One plate is positioned at h = 4cm above the other and is perfectly aligned with it. 
Both plates are clamped and initially at rest. 

The external load per unit area / is applied to the lower plate Hi only, in the upward normal 
direction. It is defined as / = 10 5 f\(x,y) f 2 (t), where 


fi{x,y) = e -100 (( x -°- 3 ) 2+ ( y - 0 - 4 ) 2 ) 

(23) 

f 2 (t ) = H(t - 10At) e - lxl ° 4 0 - 10A ‘) 2 

(24) 
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Table I. Speed-ups for online computations. 


Contact problem 

HDM size 

ROM size 

Speed-up 

Static, obstacle type 

40,000 

II 

II 

to 

o 

865 

Dynamic, obstacle type 

40,000 

p = p\ = 100 

302 

Dynamic, two-body 

20,000 

II 

II 

O 

2,229 


and H denotes the Heaviside function. 

Each plate is discretized by 100 x 100 finite elements with 1 dof per node, resulting in a 
senri-discrete HDM with a total of 20,000 dofs (10,000 dof for each plate). This HDM is 
discretized in time using the implicit second-order BDF scheme and the constant time-step 
At = 2 x 10 _4 s. 

For this problem, p and p\ are set to p = p\ =40 and Algorithm's applied with N c = 1 (no 
parametric training) to construct a ROM of size 40. Then, this ROM is applied to the solution 
for t £ [0, 0.4 s] of the same two-body dynamic contact problem as that solved using the HDM. 
Note that for At = 2 x 10 -4 s, the time-interval [0,0.4s] is sampled in 2,000 time-steps. 

Figure displays the time-histories of the HDM and ROM solutions of this problem 
at (x,y) = (0.5, 0.5). It shows that the solution delivered by the constructed ROM tracks 
remarkably well that computed using the HDM. 

Figure [12] complements Figure [III by focusing on the computed HDM and ROM solutions at 
a single time-instance, t = 0.4s, but the entire computational domain. It confirms the excellent 
accuracy of the constructed ROM. 

To illustrate the effect on the performance of a ROM of a down-sampling in time of the 
underlying HDM, Algorithm [2] is applied again to the construction of a series of contact 
ROMs in which the percentage of computed solution snapshots that are collected is uniformly 
decreased. To this effect, Figure [13] reports the variation of the relative error of the ROM 
solution with the percentage of computed solution snapshots that are collected. The reader 
can observe that overall, the relative error of the ROM solution is insensitive to a down- 
sampling in time, as long as more than 10% of the computed solution snapshots are collected 
for data compression. Beyond this limit, the relative error of the ROM solution increases 
sharply to reach the level of 52.8% when only 25 equally spaced computed solution snapshots 
are collected. 

5-4- Computational speed-up 

All model contact problems discussed above were solved in MATLAB using the quadratic 
program solver quadprog. Sparsity was accounted for in all algebraic entities of all HDMs. In 
particular, the interior-point-convex algorithm was used to solve all systems of equations 
arising from all HDMs. On the other hand, all small-scale and dense systems of equations 
arising from all ROMs were solved using the active-set algorithm. Numerical experiments 
revealed that among all algorithms available in quadprog, the aforementioned solvers are those 
which offer the best performance for the considered problems. 

All CPU times were measured using the tic-toc function on a single computational thread 
via the -singleCompThread start-up option. For each considered contact problem, the speed¬ 
up factor delivered by its ROM for the online computations is reported in Table [T| 


6 . CONCLUSIONS 

The context of this paper is set to that of the model reduction of contact problems where 
the contact conditions are enforced using Lagrange multiplier degrees of freedom (dofs). In 
this context, constructing two separate Reduced-Order Bases (ROBs), one for the primal 
displacement dofs and one for the dual Lagrange multiplier dofs, is motivated and justified 
by the positivity condition that only the dual variables must satisfy. To this effect, it is 
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shown in this paper that using the singular value decomposition method and the non-negative 
matrix factorization method to compress displacement and Lagrange multiplier snapshots, 
respectively, leads to effective primal and dual ROBs and a promising Galerkin projection 
method for the reduction of high-dimensional contact models. For parametric contact problems, 
it is also shown that the iterative greedy approach for sampling the parameter domain during 
the training of the Reduced-Order Model (ROM) can be equipped with an error indicator 
of reasonable offline computational complexity. This error indicator is based on the residual 
associated with the non-penetration condition. The computational complexity of the resulting 
iterative sampling and ROB construction procedure is dominated by the cost of a number of 
high-dimensional simulations equal to the number of sampling iterations. Specifically, for two 
parameterized, static and dynamic, model contact problems with 40,000 dofs, and one non- 
parametric two-body dynamic contact problem with 20,000 dofs, it is shown that the model 
reduction approach proposed in this paper and outlined above delivers online computational 
speedups in the range of 300 to 2,200. These are promising results which warrant the extension 
of the proposed model reduction approach to more realistic contact problems. 
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Figure 1. Illustration of a generic two-body contact problem. 



Figure 2. Parameterized obstacle g( 7 ). 




Figure 3. Convergence of Algorithm[2]for the model, parametric, static contact problem of the obstacle 
type: error bars shown in blue correspond to the ROMs constructed by randomly sampling the 

parameter space and setting p = p\ = 20 . 
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(a) 7 = ( 0 . 6 , 0 . 6 ) 



x, (:y = 0.5) 

(b) 7 = (0.330, 0.377) (Maximum ROM error) 


Figure 4. HDM and ROM solutions of the model, parametric, static contact problem of the obstacle 

type (cut view). 



(a) HDM 


(b) ROM 


Figure 5. HDM and ROM solutions for 7 = (0.330, 0.377) of the model, parametric, static contact 
problem of the obstacle type (full view). The obstacle and membrane are represented by the blue and 

red surfaces, respectively. 
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Figure 6 . Convergence of Algorithm [ 2 ] for the model, parametric, dynamic contact problem of the 

obstacle type. 



Obstacle 
-HDM 

- ROM using NNMF of dual variable 

. ROM using SVD of dual variable 

(a) (x,y) = (0.3, 0.5). 



■ Obstacle 
-HDM 

- ROM using NNMF of dual variable 

ROM using SVD of dual variable 

(b) (x,y) = (0.5, 0.5) 


Figure 7. Time-histories of the HDM and ROM solutions for 7 = (0.6, 0.6) of the model, parametric, 

dynamic contact problem of the obstacle type. 
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t 

o Obstacle 
-HDM 

-ROM using NNMF of dual variable 

. ROM using SVD of dual variable 

(a) (x,y) = (0.3,0.5) 



t 


Obstacle 
-HDM 

-ROM using NNMF of dual variable 

ROM using SVD of dual variable 

(b) (x,y) = (0.5, 0.5) 


Figure 8 . Time-histories of the HDM and ROM solutions for 7 = (0.300,0.288) of the model, 
parametric, dynamic contact problem of the obstacle type. 



(a) HDM 


(b) ROM 


Figure 9. HDM and ROM solutions at t = 0.2, 7 = (0.300, 0.288), of the model, parametric, dynamic 
contact problem of the obstacle type (full view). The obstacle and membrane are represented by the 

blue and red surfaces, respectively. 
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Figure 10. Two-body dynamic contact problem. 



Figure 11. Time-histories of the HDM and ROM solutions of the two-body dynamic contact problem 

at (x, y) = (0.5, 0.5). 



(a) HDM (b) p,k = 40 ROM 

Figure 12. Snapshots of HDM and ROM solutions of the two-body dynamic contact problem at t = 0.4. 
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Percentage of computed snapshots that are collected 


Figure 13. Effect on ROM performance of a down-sampling in time of the primal snapshots. 
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